﻿// ******************************************************************************************************
// <copyright file="Matriz.cs" company="Caroso.inc">
//     Copyright (c) Carlos Osorio All rights reserved.
// </copyright>
// ******************************************************************************************************
// The contents of this file are subject to the Mozilla Public License Version 1.1 (the "License"); 
// you may not use this file except in compliance with the License. You may obtain a copy of the License at 
// http:// www.mozilla.org/MPL/ 
// Software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTY OF 
// ANY KIND, either express or implied. See the License for the specificlanguage governing rights and 
// limitations under the License. 
// 
// 
// The Initial Developer of this version of the Original Code is hanzzoid@gmx.net
// 
// Contributor(s): (Open source contributors should list themselves and their modifications here). 
// Change Log: 
// Date           Changed By      Notes
// 22 Nov 2010    Carlos Osorio   This is an adaptation of C# Matrix library by hanzzoid@gmx.net July 2, 2007
// ******************************************************************************************************
namespace GeostatisticalTool.Classes
{
    using System;
    using System.Collections;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;
    using System.Windows.Forms;

    /// <summary>
    /// Create Matrix
    /// </summary>
    public class Matriz
    {
        #region propierties

        /// <summary>
        /// Contains the rows of the Matriz as elements, which
        /// are ArrayLists as well.
        /// </summary>
        private ArrayList values;

        /// <summary>
        /// Number of rows of the Matriz.
        /// </summary>
        private int rowCount;

        /// <summary>
        /// Number of columns of the Matriz.
        /// </summary>
        private int columnCount;

        #endregion

        #region Constructors
        /// <summary>
        /// Initializes a new instance of the Matriz class
        /// </summary>
        public Matriz()
        {
            this.values = new ArrayList();
            this.rowCount = 0;
            this.columnCount = 0;
        }
        
        /// <summary>
        /// Initializes a new instance of the Matriz class
        /// </summary>
        /// <param name="x">Complex number</param>
        public Matriz(Complex x)
        {
            this.rowCount = 1;
            this.columnCount = 1;

            this.values = new ArrayList(1);

            this.values.Add(new ArrayList(1));

            ((ArrayList)this.values[0]).Add(x);
        }
        
        /// <summary>
        /// Initializes a new instance of the Matriz class
        /// </summary>
        /// <param name="values">Array bidimension</param>
        public Matriz(double[,] values)
        {
            if (values == null)
            {
                this.values = new ArrayList();
                this.columnCount = 0;
                this.rowCount = 0;
            }

            this.rowCount = (int)values.GetLongLength(0);
            this.columnCount = (int)values.GetLongLength(1);

            this.values = new ArrayList(this.rowCount);

            for (int i = 0; i < this.rowCount; i++)
            {
                this.values.Add(new ArrayList(this.columnCount));

                for (int j = 0; j < this.columnCount; j++)
                {
                    ((ArrayList)this.values[i]).Add(new Complex(values[i, j]));
                }
            }
        }

        /// <summary>
        /// Initializes a new instance of the Matriz class
        /// </summary>
        /// <param name="values">Array simple</param>
        public Matriz(double[] values)
        {
            if (values == null)
            {
                this.values = new ArrayList();
                this.columnCount = 0;
                this.rowCount = 0;
            }

            this.rowCount = values.Length;
            this.columnCount = 1;

            this.values = new ArrayList(this.rowCount);

            for (int i = 0; i < this.rowCount; i++)
            {
                this.values.Add(new ArrayList(1));

                ((ArrayList)this.values[i]).Add(new Complex(values[i]));
            }
        }

         /// <summary>
        /// Initializes a new instance of the Matriz class  from 2-d Complex array.
        /// </summary>
        /// <param name="values"> 2-d Complex array</param>
        public Matriz(Complex[,] values)
        {
            if (values == null)
            {
                this.values = new ArrayList();
                this.columnCount = 0;
                this.rowCount = 0;
            }

            this.rowCount = (int)values.GetLongLength(0);
            this.columnCount = (int)values.GetLongLength(1);

            this.values = new ArrayList(this.rowCount);

            for (int i = 0; i < this.rowCount; i++)
            {
                this.values.Add(new ArrayList(this.columnCount));

                for (int j = 0; j < this.columnCount; j++)
                {
                    ((ArrayList)this.values[i]).Add(values[i, j]);
                }
            }
        }
       
        /// <summary>
        /// Initializes a new instance of the Matriz class
        /// </summary>
        /// <param name="n">size matrix</param>
        public Matriz(int n)
        {
            this.rowCount = n;
            this.columnCount = n;

            this.values = new ArrayList(n);

            for (int i = 0; i < n; i++)
            {
                this.values.Add(new ArrayList(n));

                for (int j = 0; j < n; j++)
                {
                    ((ArrayList)this.values[i]).Add(Complex.Zero);
                }
            }
        }

        /// <summary>
        /// Initializes a new instance of the Matriz class  by n Matriz filled with zeros; same as Zeros(m, n).
        /// </summary>
        /// <param name="m">Number of rows</param>
        /// <param name="n">Number of columns</param>
        public Matriz(int m, int n)
        {
            this.rowCount = m;
            this.columnCount = n;

            this.values = new ArrayList(m);

            for (int i = 0; i < m; i++)
            {
                this.values.Add(new ArrayList(n));

                for (int j = 0; j < n; j++)
                {
                    ((ArrayList)this.values[i]).Add(Complex.Zero);
                }
            }
        }

        #endregion

        #region fields

        /// <summary>
        /// Gets number of rows of the Matriz.
        /// </summary>
        public int RowCount
        {
            get { return this.rowCount; }
        }

        /// <summary>
        /// Gets number of columns of the Matriz.
        /// </summary>
        public int ColumnCount
        {
            get { return this.columnCount; }
        }

        #endregion

        #region indexers

        /// <summary>
        /// Access the component in row i, column j of a non-empty Matriz.
        /// </summary>
        /// <param name="i">One-based row index.</param>
        /// <param name="j">One-based column index.</param>
        /// <returns>Complex number</returns>
        public virtual Complex this[int i, int j]
        {
            get
            {
                if (i > 0 && i <= this.rowCount && j > 0 && j <= this.columnCount)
                {
                    ////Complex buf;
                    ////
                    ////try
                    ////{
                    ////    buf = (Complex)(((ArrayList)Values[i - 1])[j - 1]);
                    ////}
                    ////catch
                    ////{
                    ////    buf = new Complex((double)((int)(((ArrayList)Values[i - 1])[j - 1])));
                    ////}
                    ////
                    //// return buf;                    

                    return (Complex)((ArrayList)this.values[i - 1])[j - 1];
                }
                else
                {
                    throw new ArgumentOutOfRangeException("Indices must not exceed size of Matriz.");
                }
            }

            set
            {
                if (i <= 0 || j <= 0)
                {
                    throw new ArgumentOutOfRangeException("Indices must be real positive.");
                }

                if (i > this.rowCount)
                {
                    // dynamically add i-Rows new rows...
                    for (int k = 0; k < i - this.rowCount; k++)
                    {
                        this.values.Add(new ArrayList(this.columnCount));

                        // ...with Cols columns
                        for (int t = 0; t < this.columnCount; t++)
                        {
                            ((ArrayList)this.values[this.rowCount + k]).Add(Complex.Zero);
                        }
                    }

                    this.rowCount = i; // ha!
                }

                if (j > this.columnCount)
                {
                    //// dynamically add j-Cols columns to each row
                    for (int k = 0; k < this.rowCount; k++)
                    {
                        for (int t = 0; t < j - this.columnCount; t++)
                        {
                            ((ArrayList)this.values[k]).Add(Complex.Zero);
                        }
                    }

                    this.columnCount = j;
                }

                ((ArrayList)this.values[i - 1])[j - 1] = value;
                ////this.Values[i - 1, j - 1] = value; 
            }
        }

        /// <summary>
        /// Access to the i-th component of an n by one Matriz (column vector)
        /// or one by n Matriz (row vector).
        /// </summary>
        /// <param name="i">One-based index.</param>
        /// <returns>Complex number</returns>
        public virtual Complex this[int i]
        {
            get
            {
                ////Complex buf;

                ////if (this.rowCount == 1)
                ////    try
                ////    {
                ////        buf = (Complex)(((ArrayList)Values[0])[i - 1]);
                ////    }
                ////    catch
                ////    {
                ////        buf = new Complex((double)((int)(((ArrayList)Values[0])[i - 1])));
                ////    }
                ////else
                ////    try
                ////    {
                ////        buf = (Complex)(((ArrayList)Values[i - 1])[0]);
                ////    }
                ////    catch
                ////    {
                ////        buf = new Complex((double)((int)(((ArrayList)Values[i - 1])[0])));
                ////    }

                ////return buf;
                //// row vector
                if (this.RowCount == 1) 
                {
                    return (Complex)((ArrayList)this.values[0])[i - 1];
                }
                else if (this.ColumnCount == 1) 
                {
                    //// coumn vector
                    return (Complex)((ArrayList)this.values[i - 1])[0];
                }
                else 
                {
                    //// neither
                    throw new InvalidOperationException("General Matriz acces requires double indexing.");
                }
            }

            set
            {
                if (this.rowCount == 1)
                {
                    // row vector

                    // dynamically extend vector if necessary
                    if (i > this.columnCount)
                    {
                        // dynamically add j-Cols columns to each row
                        for (int t = 0; t < i - this.columnCount; t++)
                        {
                            ((ArrayList)this.values[0]).Add(Complex.Zero);
                        }

                        this.columnCount = i;
                    }

                    ((ArrayList)this.values[0])[i - 1] = value;
                }
                else if (this.columnCount == 1)
                {
                    //// column vector

                    if (i > this.rowCount)
                    {
                        // dynamically add i-Rows new rows...
                        for (int k = 0; k < i - this.rowCount; k++)
                        {
                            this.values.Add(new ArrayList(this.columnCount));

                            //// ...with one column each
                            ((ArrayList)this.values[this.rowCount + k]).Add(Complex.Zero);
                        }

                        this.rowCount = i;
                        //// ha!
                    }

                    ((ArrayList)this.values[i - 1])[0] = value;
                }
                else
                {
                    throw new InvalidOperationException("Cannot access multidimensional Matriz via single index.");
                }
            }
        }

        #endregion

        #region static methods

        /// <summary>
        /// Generates diagonal Matriz
        /// </summary>
        /// <param name="diag_vector">column vector containing the diag elements</param>
        /// <returns>Diagonal matrix</returns>
        public static Matriz Diag(Matriz diag_vector)
        {
            int dim = diag_vector.VectorLength();

            if (dim == 0)
            {
                throw new ArgumentException("diag_vector must be 1xN or Nx1");
            }

            Matriz m = new Matriz(dim, dim);

            for (int i = 1; i <= dim; i++)
            {
                m[i, i] = diag_vector[i];
            }

            return m;
        }

        /// <summary>
        /// Generates diagonal Matriz
        /// </summary>
        /// <param name="diag_vector">column vector containing the diag elements</param>
        /// <param name="offset">Offset of matrix</param>
        /// <returns>Diagonal matrix</returns>
        public static Matriz Diag(Matriz diag_vector, int offset)
        {
            int dim = diag_vector.VectorLength();

            if (dim == 0)
            {
                throw new ArgumentException("diag_vector must be 1xN or Nx1.");
            }
            ////if (Math.Abs(offset) >= dim)
            ////    throw new ArgumentException("Absolute value of offset must be less than length of diag_vector.");

            Matriz m = new Matriz(dim + Math.Abs(offset), dim + Math.Abs(offset));
            dim = m.RowCount;

            if (offset >= 0)
            {
                for (int i = 1; i <= dim - offset; i++)
                {
                    m[i, i + offset] = diag_vector[i];
                }
            }
            else
            {
                for (int i = 1; i <= dim + offset; i++)
                {
                    m[i - offset, i] = diag_vector[i];
                }
            }

            return m;
        }

        /// <summary>
        /// Implements the dot product of two vectors.
        /// </summary>
        /// <param name="v">Row or column vector.</param>
        /// <param name="w">Row or column of other vector.</param>
        /// <returns>Dot product.</returns>
        public static Complex Dot(Matriz v, Matriz w)
        {
            int m = v.VectorLength();
            int n = w.VectorLength();

            if (m == 0 || n == 0)
            {
                throw new ArgumentException("Arguments need to be vectors.");
            }
            else if (m != n)
            {
                throw new ArgumentException("Vectors must be of the same length.");
            }

            Complex buf = Complex.Zero;

            for (int i = 1; i <= m; i++)
            {
                buf += v[i] * w[i];
            }

            return buf;
        }
       
        /// <summary>
        /// Creates n by n identity Matriz.
        /// </summary>
        /// <param name="n">Number of rows and columns respectively.</param>
        /// <returns>n by n identity Matriz.</returns>
        public static Matriz Identity(int n)
        {
            return Diag(Ones(n, 1));
        }

        /// <summary>
        /// Creates m by n Matriz filled with ones.
        /// </summary>
        /// <param name="m">Number of rows.</param>
        /// <param name="n">Number of columns.</param>
        /// <returns>m by n Matriz filled with ones.</returns>        
        public static Matriz Ones(int m, int n)
        {
            Matriz mat = new Matriz(m, n);

            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    ((ArrayList)mat.values[i])[j] = Complex.One;
                }
            }

            return mat;
        }

        /// <summary>
        /// Creates n by n Matriz filled with ones.
        /// </summary>        
        /// <param name="n">Number of columns.</param>
        /// <returns>n by n Matriz filled with ones.</returns>        
        public static Matriz Ones(int n)
        {
            Matriz m = new Matriz(n);

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    ((ArrayList)m.values[i])[j] = Complex.One;
                }
            }

            return m;
        }

        /// <summary>
        /// Define operator equals
        /// </summary>
        /// <param name="a">First matrix to be compared</param>
        /// <param name="b">Second matrix to be compared</param>
        /// <returns>True or false</returns>
        public static bool operator ==(Matriz a, Matriz b)
        {
            if (a.RowCount != b.RowCount || a.ColumnCount != b.ColumnCount)
            {
                return false;
            }

            for (int i = 1; i <= a.RowCount; i++)
            {
                for (int j = 1; j <= a.ColumnCount; j++)
                {
                    if (a[i, j] != b[i, j])
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        /// <summary>
        /// Define operator differt
        /// </summary>
        /// <param name="a">First matrix to be compared</param>
        /// <param name="b">Second matrix to be compared</param>
        /// <returns>True or false</returns>
        public static bool operator !=(Matriz a, Matriz b)
        {
            return !(a == b);
        }

        /// <summary>
        /// Define operator addition
        /// </summary>
        /// <param name="a">First matrix in the operation</param>
        /// <param name="b">Second matrix in the operation</param>
        /// <returns>Matrix result</returns>
        public static Matriz operator +(Matriz a, Matriz b)
        {
            if (a.RowCount != b.RowCount || a.ColumnCount != b.ColumnCount)
            {
                throw new ArgumentException("Matrices must be of the same dimension.");
            }

            for (int i = 1; i <= a.RowCount; i++)
            {
                for (int j = 1; j <= a.ColumnCount; j++)
                {
                    a[i, j] += b[i, j];
                }
            }

            return a;
        }

        /// <summary>
        /// Define operator substraction
        /// </summary>
        /// <param name="a">First matrix in the operation</param>
        /// <param name="b">Second matrix in the operation</param>
        /// <returns>Matrix result</returns>
        public static Matriz operator -(Matriz a, Matriz b)
        {
            if (a.RowCount != b.RowCount || a.ColumnCount != b.ColumnCount)
            {
                throw new ArgumentException("Matrices must be of the same dimension.");
            }

            for (int i = 1; i <= a.RowCount; i++)
            {
                for (int j = 1; j <= a.ColumnCount; j++)
                {
                    a[i, j] -= b[i, j];
                }
            }

            return a;
        }

        /// <summary>
        /// Define operator substraction
        /// </summary>
        /// <param name="a">First matrix in the operation</param>
        /// <returns>Matrix result</returns>
        public static Matriz operator -(Matriz a)
        {
            for (int i = 1; i <= a.RowCount; i++)
            {
                for (int j = 1; j <= a.ColumnCount; j++)
                {
                    a[i, j] = -a[i, j];
                }
            }

            return a;
        }

        /// <summary>
        /// Define operator multiplication
        /// </summary>
        /// <param name="a">First matrix in the operation</param>
        /// <param name="b">Second matrix in the operation</param>
        /// <returns>Matrix result</returns>
        public static Matriz operator *(Matriz a, Matriz b)
        {
            if (a.ColumnCount != b.RowCount)
            {
                throw new ArgumentException("Inner Matriz dimensions must agree.");
            }

            Matriz c = new Matriz(a.RowCount, b.ColumnCount);

            for (int i = 1; i <= a.RowCount; i++)
            {
                for (int j = 1; j <= b.ColumnCount; j++)
                {
                    c[i, j] = Dot(a.Row(i), b.Column(j));
                }
            }

            return c;
        }

        /// <summary>
        /// Define operator multiplication with a complex number
        /// </summary>
        /// <param name="a">First matrix in the operation</param>
        /// <param name="x">Complex number in the operation</param>
        /// <returns>Matrix result</returns>
        public static Matriz operator *(Matriz a, Complex x)
        {
            Matriz b = new Matriz(a.rowCount, a.columnCount);

            for (int i = 1; i <= a.RowCount; i++)
            {
                for (int j = 1; j <= a.ColumnCount; j++)
                {
                    b[i, j] = a[i, j] * x;
                }
            }

            return b;
        }

        /// <summary>
        /// Define operator multiplication with a complex number 
        /// </summary>
        /// <param name="x">First complex number in the operation</param>
        /// <param name="a">Second matrix in the operation</param>
        /// <returns>Matrix result</returns>
        public static Matriz operator *(Complex x, Matriz a)
        {
            Matriz b = new Matriz(a.RowCount, a.ColumnCount);

            for (int i = 1; i <= a.RowCount; i++)
            {
                for (int j = 1; j <= a.ColumnCount; j++)
                {
                    b[i, j] = a[i, j] * x;
                }
            }

            return b;
        }

        /// <summary>
        /// Define operator multiplication with a double number 
        /// </summary>
        /// <param name="a">First Matrix in the operation</param>
        /// <param name="x">Second double number </param>
        /// <returns>Matrix result</returns>
        public static Matriz operator *(Matriz a, double x)
        {
            return (new Complex(x)) * a;
        }

        /// <summary>
        /// Define operator multiplication with a double number 
        /// </summary>
        /// <param name="x">First double number</param>
        /// <param name="a">Matrix in the operation</param>
        /// <returns>Matrix result</returns>
        public static Matriz operator *(double x, Matriz a)
        {
            return (new Complex(x)) * a;
        }

        /// <summary>
        /// Define operator division with a complex number 
        /// </summary>
        /// <param name="a">First Matrix in the operation</param>
        /// <param name="x">Second double number </param>
        /// <returns>Matrix result</returns>
        public static Matriz operator /(Matriz a, Complex x)
        {
            return (1 / x) * a;
        }

        /// <summary>
        /// Define operator division with a double number 
        /// </summary>
        /// <param name="a">First double number</param>
        /// <param name="x">Second Matrix in the operation </param>
        /// <returns>Matrix result</returns>
        public static Matriz operator /(Matriz a, double x)
        {
            return (new Complex(1 / x)) * a;
        }

        /// <summary>
        /// Define operator power 
        /// </summary>
        /// <param name="a">First Matrix in the operation</param>
        /// <param name="k">Second integer number </param>
        /// <returns>Matrix result</returns>
        public static Matriz operator ^(Matriz a, int k)
        {
            if (k < 0)
            {
                if (a.IsSquare())
                {
                    return a.InverseLeverrier() ^ (-k);
                }
                else
                {
                    throw new InvalidOperationException("Cannot take non-square Matriz to the power of zero.");
                }
            }
            else if (k == 0)
            {
                if (a.IsSquare())
                {
                    return Matriz.Identity(a.RowCount);
                }
                else
                {
                    throw new InvalidOperationException("Cannot take non-square Matriz to the power of zero.");
                }
            }
            else if (k == 1)
            {
                if (a.IsSquare())
                {
                    return a;
                }
                else
                {
                    throw new InvalidOperationException("Cannot take non-square Matriz to the power of one.");
                }
            }
            else
            {
                Matriz m = a;
                for (int i = 1; i < k; i++)
                {
                    m *= a;
                }

                return m;
            }
        }

        /// <summary>
        /// Creates m by n Matriz filled with zeros.
        /// </summary>
        /// <param name="m">Number of rows.</param>
        /// <param name="n">Number of columns.</param>
        /// <returns>m by n Matriz filled with zeros.</returns>
        public static Matriz Zeros(int m, int n)
        {
            return new Matriz(m, n);
        }
        
        #endregion
        
        #region methods

        /// <summary>
        /// Provides a shallow copy of this Matriz in O(m).
        /// </summary>
        /// <returns>Matrix clonned</returns>
        public Matriz Clone()
        {
            Matriz a = new Matriz();
            a.rowCount = this.rowCount;
            a.columnCount = this.columnCount;

            for (int i = 0; i < this.rowCount; i++)
            {
                a.values.Add(((ArrayList)this.values[i]).Clone());
            }

            return a;
        }

        /// <summary>
        /// Retrieves column with one-based index j.
        /// </summary>
        /// <param name="j">Index j of the matrix</param>
        /// <returns>j-th column of the matrix</returns>
        public Matriz Column(int j)
        {
            Matriz buf = new Matriz(this.rowCount, 1);

            for (int i = 1; i <= this.rowCount; i++)
            {
                buf[i] = this[i, j];
            }

            return buf;
        }

        /// <summary>
        /// Conjuagtes and transposes a Matriz.
        /// </summary>
        /// <returns>Matrix transpose</returns>
        public Matriz ConjTranspose()
        {
            return this.Transpose().Conjugate();
        }

        /// <summary>
        /// Replaces each Matriz entry z = x + iy with x - iy.
        /// </summary>
        /// <returns>Conjugated Matriz.</returns>
        public Matriz Conjugate()
        {
            Matriz m = new Matriz(this.rowCount, this.columnCount);

            for (int i = 1; i <= this.rowCount; i++)
            {
                for (int j = 1; j <= this.columnCount; j++)
                {
                    m[i, j] = Complex.Conj(this[i, j]);
                }
            }

            return m;
        }

        /// <summary>
        /// Deletes row at specifies index.
        /// </summary>
        /// <param name="i">One-based index at which to delete.</param>
        public void DeleteRow(int i)
        {
            if (i <= 0 || i > this.rowCount)
            {
                throw new ArgumentException("Index must be positive and <= number of rows.");
            }

            this.values.RemoveAt(i - 1);
            this.rowCount--;
        }

        /// <summary>
        /// Deletes column at specifies index.
        /// </summary>
        /// <param name="j">One-based index at which to delete.</param>
        public void DeleteColumn(int j)
        {
            if (j <= 0 || j > this.columnCount)
            {
                throw new ArgumentException("Index must be positive and <= number of cols.");
            }

            for (int i = 0; i < this.rowCount; i++)
            {
                ((ArrayList)this.values[i]).RemoveAt(j - 1);
            }

            this.columnCount--;
        }

        /// <summary>
        /// Calcs determinant of square Matriz
        /// </summary>
        /// <returns>Complex number</returns>
        public Complex Determinant()
        {
            if (!this.IsSquare())
            {
                throw new InvalidOperationException("Cannot calc determinant of non-square Matriz.");
            }
            else if (this.columnCount == 1)
            {
                return this[1, 1];
            }
            else if (this.IsTrapeze()) 
            {
                //// is square, therefore triangular
                return this.DiagProd();
            }
            else
            {
                //// perform LU-decomposition & return product of diagonal elements of U
                Matriz x = this.Clone();

                //// for speed concerns, use this
                //x..LU();
                //return x.DiagProd();

                //// this is slower and needs more memory... .
                Matriz p = x.LUSafe();
                return (double)p.Signum() * x.DiagProd();
            }
        }

        /// <summary>
        /// Computes product of main diagonal entries.
        /// </summary>
        /// <returns>Product of diagonal elements</returns>
        public Complex DiagProd()
        {
            Complex buf = Complex.One;
            int dim = Math.Min(this.rowCount, this.columnCount);

            for (int i = 1; i <= dim; i++)
            {
                buf *= this[i, i];
            }

            return buf;
        }

        /// <summary>
        /// Extracts main diagonal vector of the Matriz as a column vector.
        /// </summary>
        /// <returns>Vector diagonal</returns>
        public Matriz DiagVector()
        {
            if (!this.IsSquare())
            {
                throw new InvalidOperationException("Cannot get diagonal of non-square Matriz.");
            }

            Matriz v = new Matriz(this.columnCount, 1);

            for (int i = 1; i <= this.columnCount; i++)
            {
                v[i] = this[i, i];
            }

            return v;
        }

        /// <summary>
        /// Override Equal comparison
        /// </summary>
        /// <param name="obj">Object to be compared</param>
        /// <returns>True or False</returns>
        public override bool Equals(object obj)
        {
            return obj.ToString() == this.ToString();
        }

        /// <summary>
        /// Inverts square Matriz as long as det != 0.
        /// </summary>
        /// <returns>Inverse of Matriz.</returns>
        public Matriz Inverse()
        {
            if (!this.IsSquare())
            {
                throw new InvalidOperationException("Cannot invert non-square Matriz.");
            }

            Complex det = this.Determinant();

            if (det == Complex.Zero)
            {
                throw new InvalidOperationException("Cannot invert (nearly) singular Matriz.");
            }

            int n = this.columnCount;

            if (n == 1)
            {
                return new Matriz(1 / det);
            }

            if (this.IsReal() && this.IsOrthogonal())
            {
                return this.Transpose();
            }
            else if (this.IsUnitary())
            {
                return this.ConjTranspose();
            }

            if (this.IsDiagonal())
            {
                Matriz d = this.DiagVector();

                for (int i = 1; i <= n; i++)
                {
                    d[i] = 1 / d[i];
                }

                return Diag(d);
            }

            Complex[,] buf = new Complex[n, n];

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    buf[i, j] = Math.Pow(-1, i + j) * this.Minor(j + 1, i + 1).Determinant();
                }
            }

            return new Matriz(buf) / det;
        }

         /// <summary>
        /// Alternative matrix inversion using Leverrier's formula
        /// </summary>
        /// <returns>Inverse of matrix.</returns>
        public Matriz InverseLeverrier()
        {
            if (!this.IsSquare())
            {
                throw new InvalidOperationException("Cannot invert non-square matrix.");
            }
            ////else if (this.Determinant() == 0)
            ////    throw new InvalidOperationException("Cannot invert (nearly) singular matrix.");

            int n = this.rowCount;
            Matriz id = Identity(n);
            Matriz b = id;
            Complex alpha;

            for (int k = 1; k < n; k++)
            {                
                Matriz buf = this * b; // DEBUG                
                Complex buf2 = buf.Trace(); // DEBUG
                alpha = ((double)1 / k) * buf.Trace();
                b = (alpha * id) - buf;
            }

            Matriz buf3 = this * b; // DEBUG                
            Complex buf4 = buf3.Trace(); // DEBUG
            alpha = (this * b).Trace() / n;
            if (alpha != Complex.Zero)
            {
                return b / alpha;
            }
            else
            {
                throw new InvalidOperationException("WARNING: Matrix nearly singular or badly scaled.");
            }
        }

        /// <summary>
        /// Checks if Matriz is diagonal Matriz.
        /// </summary>
        /// <returns>True iff Matriz is diagonal.</returns>
        public bool IsDiagonal()
        {
            return this.Clone() - Diag(this.DiagVector()) == Zeros(this.rowCount, this.columnCount);
        }

        /// <summary>
        /// Checks if number of rows equals number of columns.
        /// </summary>
        /// <returns>True iff Matriz is n by n.</returns>
        public bool IsSquare()
        {
            return this.columnCount == this.rowCount;
        }

        /// <summary>
        /// Checks if Matriz consists only of real entries.
        /// </summary>
        /// <returns>True iff all entries are real.</returns>
        public bool IsReal()
        {
            for (int i = 1; i <= this.rowCount; i++)
            {
                for (int j = 1; j <= this.columnCount; j++)
                {
                    if (!this[i, j].IsReal())
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        /// <summary>
        /// Checks if Matriz is lower or upper trapeze.
        /// </summary>
        /// <returns>True iff Matriz is trapeze.</returns>
        public bool IsTrapeze()
        {
            return this.IsUpperTrapeze() || this.IsLowerTrapeze();
        }

        /// <summary>
        /// Checks if Matriz is trapeze and square.
        /// </summary>
        /// <returns>True iff Matriz is triangular.</returns>
        public bool IsTriangular()
        {
            return this.IsLowerTriangular() || this.IsUpperTriangular();
        }
        
        /// <summary>
        /// Matriz is unitary, iff A^H*A = id, where A^H is the conjugated transpose of A.
        /// </summary>
        /// <returns>True iff Matriz is unitary.</returns>
        public bool IsUnitary()
        {
            if (!this.IsSquare())
            {
                return false;
            }

            return this.ConjTranspose() * this == Identity(this.rowCount);
        }

        /// <summary>
        /// Checks if Matriz is square and upper trapeze.
        /// </summary>
        /// <returns>True iff Matriz is upper triangular.</returns>
        public bool IsUpperTriangular()
        {
            return this.IsSquare() && this.IsUpperTrapeze();
        }

        /// <summary>
        /// Checks if Matriz is square and lower trapeze.
        /// </summary>
        /// <returns>True iff Matriz is lower triangular.</returns>
        public bool IsLowerTriangular()
        {
            return this.IsSquare() && this.IsLowerTrapeze();
        }
        
        /// <summary>
        /// Checks for orthogonality by testing if A*A' == id.
        /// </summary>
        /// <returns>True iff Matriz is orthogonal.</returns>
        public bool IsOrthogonal()
        {
            return this.IsSquare() && this * this.Transpose() == Identity(this.rowCount);
        }
        
        /// <summary>
        /// Checks if A in i j are equals 0 for i less than j.
        /// </summary>
        /// <returns>True iff Matriz is upper trapeze.</returns>
        public bool IsUpperTrapeze()
        {
            for (int j = 1; j <= this.columnCount; j++)
            {
                for (int i = j + 1; i <= this.rowCount; i++)
                {
                    if (this[i, j] != 0)
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        /// <summary>
        /// Checks if A[i, j] == 0 for i > j.
        /// </summary>
        /// <returns>True iff Matriz is lower trapeze.</returns>
        public bool IsLowerTrapeze()
        {
            for (int i = 1; i <= this.rowCount; i++)
            {
                for (int j = i + 1; j <= this.columnCount; j++)
                {
                    if (this[i, j] != 0)
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        /// <summary>
        /// Performs safe LU-decomposition of this instance with column pivoting 
        /// and saves L and U
        /// within, where the diagonal elements belong to U
        /// (the ones of L are ones...)
        /// </summary>
        /// <returns>Permutation Matriz P with P*this = L*U</returns>
        /// <remarks>This needs additional time O(n^2).</remarks>
        public Matriz LUSafe()
        {
            if (!this.IsSquare())
            {
                throw new InvalidOperationException("Cannot perform LU-decomposition of non-square Matriz.");
            }

            int n = this.columnCount;

            Matriz p = Identity(n); // permutation Matriz
            int m;

            for (int j = 1; j <= n; j++)
            {
                ////--> this test means probably deceleration
                ////if (j < n && this.Extract(j + 1, n, j, j) == Zeros(n - j, 1))
                ////    continue;
                //// find index m with |this[m,j]| >= |this[i,j]| for all i in {j,...,n}
                if (j < n)
                {
                    m = j;

                    for (int i = j + 1; i <= n; i++)
                    {
                        if (Complex.Abs(this[i, j]) > Complex.Abs(this[m, j]))
                        {
                            m = i;
                        }
                    }

                    if (m > j) 
                    {
                        //// <=> j2 != j
                        p.SwapRows(j, m);
                        this.SwapRows(j, m);
                    }

                    if (this[j, j] == 0)
                    {
                        throw new DivideByZeroException("Warning: Matriz close to singular.");
                    }
                }

                for (int k = 1; k < j; k++)
                {
                    for (int i = k + 1; i <= n; i++)
                    {
                        this[i, j] = this[i, j] - (this[i, k] * this[k, j]);
                    }
                }

                for (int i = j + 1; i <= n; i++)
                {
                    this[i, j] = this[i, j] / this[j, j];
                }
            }

            return p;
        }

        /// <summary>
        /// Calcs the matrix that results in the clearing of a
        /// specified row and a specified column
        /// </summary>        
        /// <param name="i">Position  column</param>
        /// <param name="j">Possition row</param>
        /// <returns>The Minor Matrix</returns>
        public Matriz Minor(int i, int j)
        {
            //// THIS IS THE LOW-LEVEL SOLUTION ~ O(n^2)
            ////Complex[,] buf = new Complex[RowCount - 1, ColumnCount - 1];
            ////int r = 0;
            ////int c = 0;

            ////for (int i = 1; i <= RowCount; i++)
            ////{
            ////    if (i != row)
            ////    {
            ////        for (int j = 1; j <= ColumnCount; j++)
            ////        {
            ////            if (j != col)
            ////            {
            ////                buf[r, c] = this[i, j];
            ////                c++;
            ////            }
            ////        }

            ////        c = 0;
            ////        r++;
            ////    }
            ////}

            ////return new Matrix(buf);  

            //// THIS IS THE HIGH-LEVEL SOLUTION ~ O(n)

            Matriz a = this.Clone();

            a.DeleteRow(i);
            a.DeleteColumn(j);

            return a;
        }

        /// <summary>
        /// Retrieves row with one-based index i.
        /// </summary>
        /// <param name="i">Define position searching</param>
        /// <returns>i-th row in the matrix</returns>
        public Matriz Row(int i)
        {
            if (i <= 0 || i > this.rowCount)
            {
                throw new ArgumentException("Index exceed matrix dimension.");
            }

            ////return (new Matrix((Complex[])((ArrayList)Values[i - 1]).ToArray(typeof(Complex)))).Transpose();

            Matriz buf = new Matriz(this.columnCount, 1);

            for (int j = 1; j <= this.columnCount; j++)
            {
                buf[j] = this[i, j];
            }

            return buf;
        }

        /// <summary>
        /// Computes signum of a permutation Matriz, which is 1 for an even
        /// number of swaps and -1 for an odd number of swaps. WARNING: 
        /// if *this is not a permutation Matriz (e.i. a permutation of Id),
        /// garbage is returned.
        /// </summary>
        /// <returns>double value</returns>
        public double Signum()
        {
            double buf = 1;

            int n = this.rowCount;
            double fi, fj;

            for (int i = 1; i < n; i++)
            {
                for (fi = 1; fi < n && this[i, (int)fi] != Complex.One; fi++)
                {
                }

                for (int j = i + 1; j <= n; j++)
                {
                    for (fj = 1; fj <= n && this[j, (int)fj] != Complex.One; fj++)
                    {
                    }

                    buf *= (fi - fj) / (i - j);
                }
            }

            return buf;
        }

        /// <summary>
        /// Swaps rows at specified indices. The latter do not have to be ordered.
        /// When equal, nothing is done.
        /// </summary>
        /// <param name="i1">One-based index of first row.</param>
        /// <param name="i2">One-based index of second row.</param>        
        public void SwapRows(int i1, int i2)
        {
            if (i1 <= 0 || i1 > this.rowCount || i2 <= 0 || i2 > this.rowCount)
            {
                throw new ArgumentException("Indices must be positive and <= number of rows.");
            }

            if (i1 == i2)
            {
                return;
            }

            ArrayList buf = (ArrayList)this.values[--i1];
            this.values[i1] = this.values[--i2];
            this.values[i2] = buf;
        }
        
        /// <summary>
        /// Redefine the funtion GetHashCode
        /// </summary>
        /// <returns>Integer of the Haschode</returns>
        public override int GetHashCode()
        {
            return -1;
        }

        /// <summary>
        /// Override the function ToString
        /// </summary>
        /// <returns>String of the matrix</returns>
        public override string ToString()
        {
            string s = string.Empty;
            Complex buf;
            Complex aux;
            aux = new Complex(double.NaN, double.NaN);
           //// MessageBox.Show(aux.ToString());

            for (int i = 1; i <= this.rowCount; i++)
            {
                for (int j = 1; j <= this.columnCount; j++)
                {
                    buf = this[i, j];
                    
                    if (buf.Realnumber == double.PositiveInfinity || buf.Realnumber == double.NegativeInfinity || buf.Imaginarynumber == double.PositiveInfinity || buf.Imaginarynumber == double.NegativeInfinity)
                    {
                        s += "oo";
                    }
                    else if (buf.Realnumber == double.NaN || buf.Imaginarynumber == double.NaN)
                    {
                        s += "?";
                    }
                    else
                    {
                        s += buf.ToString();
                    }

                    s += ";" + "\t";
                }

                s += "\\" + System.Environment.NewLine;
            }

            return s;
        }

        /// <summary>
        /// Export to array
        /// </summary>
        /// <returns>External array</returns>
        public double[,] ToArray()
        {
            double[,] mat = new double[this.rowCount, this.rowCount];
            Complex buf;
            Complex aux;
            aux = new Complex(double.NaN, double.NaN);

            for (int i = 1; i <= this.rowCount; i++)
            {
                for (int j = 1; j <= this.columnCount; j++)
                {
                    buf = this[i, j];
                    if (buf.Equals(aux))
                    {
                        mat[i - 1, j - 1] = 0;
                    }
                    else
                    {
                        mat[i - 1, j - 1] = buf.Realnumber;
                    }
                }
            }

            return mat;
        }

        /// <summary>
        /// Calcs trace of the matrix.
        /// </summary>
        /// <returns>Sum of diagonal elements.</returns>
        public Complex Trace()
        {
            if (!this.IsSquare())
            {
                throw new InvalidOperationException("Cannot calc trace of non-square matrix.");
            }

            Complex buf = Complex.Zero;

            for (int i = 1; i <= this.rowCount; i++)
            {
                buf += this[i, i];
            }

            return buf;
        }

        /// <summary>
        /// Swaps each Matriz entry A[i, j] with A[j, i].
        /// </summary>
        /// <returns>A transposed Matriz.</returns>
        public Matriz Transpose()
        {
            Matriz m = new Matriz(this.columnCount, this.rowCount);

            for (int i = 1; i <= this.columnCount; i++)
            {
                for (int j = 1; j <= this.rowCount; j++)
                {
                    m[i, j] = this[j, i];
                }
            }

            return m;
        }
        
        /// <summary>
        /// Checks if Matriz is n by one or one by n.
        /// </summary>
        /// <returns>Length, if vector; zero else.</returns>
        public int VectorLength()
        {
            if (this.columnCount > 1 && this.rowCount > 1)
            {
                return 0;
            }
            else
            {
                return Math.Max(this.columnCount, this.rowCount);
            }
        }

        #endregion
    }
}
